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Simulations are carried out based on the dynamical mean-field theory (DMFT) in order to investi¬ 
gate the properties of correlated thin films for various values of the chemical potential, temperature, 
interaction strength, and applied transverse electric field. Application of a sufficiently strong field 
to a thin film at half-filling leads to the appearance of conducting regions near the surfaces of the 
film, whereas in doped slabs the application of a field leads to a conductivity enhancement on one 
side of the film and a gradual transition to the insulating state on the opposite side. In addition 
to the inhomogeneous DMFT, an independent layer approximation (ILA) is considered, in which 
the properties of each layer are approximated by a homogeneous bulk environment. A comparison 
between the two approaches reveals that the less expensive ILA results are in good agreement with 
the DMFT approach, except close to the metal-to-insulator transition points and in the layers im¬ 
mediately at the film surfaces. The hysteretic behavior (memory effect) characteristic of the bulk 
doping driven Mott transition persists in the slab. 


I. INTRODUCTION 

Due to their interesting and useful properties, strongly 
correlated materials have attracted the attention of both 
theorists 1-3 and experimentalists 4-8 in recent decades. In 
particular, the discoveries of high-temperature supercon¬ 
ductivity in cuprates and of colossal magnetoresistance 
in manganites spurred interest in the class of materi¬ 
als exhibiting the Mott metal-to-insulator transition. 9,10 
The existence of this correlation-driven transition in bulk 
transition-metal oxide systems has encouraged investiga¬ 
tions for potential applications in electronics. 5- ' Intrigu¬ 
ing physics has also been observed at interfaces between 
strongly correlated materials; 11 for example, a metallic 
or even superconducting phase appears at the interface 
between LaTiC >3 and SrTiOa, which are insulating para- 
magnets in bulk. 8,12 The appearance of such interface 
phases can be partially understood at a qualitative level 
as the result of charge transfer. A detailed understand¬ 
ing, however, remains elusive. Due to the difficulties in¬ 
herent in the theoretical treatment of strongly-correlated 
systems, attaining a detailed understanding may be quite 
challenging even when dealing with simpler bulk systems. 

The exponential growth of the Hilbert space with 
the number of particles makes a direct numerical solu¬ 
tion of an interacting quantum system unfeasible, ex¬ 
cept for systems with a very limited number of parti¬ 
cles. Standard theoretical approaches such as density 
functional theory (DFT) are known to fail when electron- 
electron interactions are strong. 13 This has stimulated 
the development of alternative methods. The Gutzwiller 
approximation 14-16 and slave boson techniques 17,18 have 
been used to treat the low-energy physics of bulk materi¬ 
als and - more recently - thin films 19-24 in the presence of 
an electric field. 21,22 The slave boson approach has been 
extended to allow access to higher-energy excitations 25 , 
while the Gutzwiller approximation only provides insight 


into the metallic state. Other notable approaches include 
the density renormalization group (DMRG), which works 
best for one-dimensional systems. 26 Over the last 25 years 
dynamical mean-field theory (DMFT) has emerged as 
one of the most promising frameworks for the treat¬ 
ment of strongly correlated systems. 2 '’ 28 DMFT involves 
a mapping of a lattice problem, such as the Hubbard 
model, 29 onto an interacting quantum impurity model 30 
supplemented by a self-consistency condition. The map¬ 
ping is exact in the limit of infinite dimensions. An 
advantage of DMFT is that it is formulated in the 
thermodynamic limit and that it allows one to non- 
perturbatively interpolate between the strong and weak 
coupling regimes, thereby treating different phases on an 
equal footing. In its simplest version, single-site DMFT, 
it neglects spatial correlations, but captures local tem¬ 
poral correlations. In combination with the local density 
approximation (LDA) to DFT it has provided valuable 
new insights into the physics of a number of strongly cor¬ 
related bulk materials. 13 ’ 31,32 

A better understanding of the metal-insulator transi¬ 
tion in systems with lower translational symmetry (such 
as thin films or interfaces) and in the presence of an 
electric field is likely to aid experimental efforts to con¬ 
trol the properties of such systems, which is relevant for 
possible applications in e.g. memory devices and high¬ 
speed electronics. Potthoff and Nolting 33,34 used a spa¬ 
tially resolved version of DMFT (inhomogeneous DMFT, 
or IDMFT) to study layered systems. Similar tech¬ 
niques were later applied to inhomogeneous multilayered 
nanostructures. 1,2,35,36 DMFT-based approaches were re¬ 
cently proposed for the study of "Mott p-n junctions" 37 
and correlated capacitors. 38 In these approaches each 
layer was approximated by a two-dimensional bulk sys¬ 
tem with the appropriate value of the chemical potential, 
and a Poisson solver was used to study the charge redis¬ 
tribution in different regimes. 

In this paper we use the IDMFT to investigate the ef- 
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feet of an applied electric field on the properties of a slab 
(a stack of coupled correlated two-dimensional layers) at 
various temperatures and doping levels. We also approx¬ 
imate the properties of each layer in the slab by a three- 
dimensional bulk system in which the value of the chem¬ 
ical potential is equal to the local value of the potential 
in that layer. We refer to this approach as the indepen¬ 
dent layer approach (ILA) and we carry out a systematic 
comparison between the ILA and the full DMFT calcula¬ 
tion. This comparison is of some interest, since the ILA 
is significantly faster than the full calculation. We inves¬ 
tigate the effect of screening by considering the Coulomb 
interaction between planes in the slab at a mean-field 
level. 

The structure of this paper is as follows. We describe 
the model and method in Sec. II. Our results are dis¬ 
cussed in Sec. III. Sec. IV provides a summary. 


II. MODEL AND METHOD 

The general approach that we use is described in detail 
elsewhere , 33 ’ 34 and we only give a brief description here. 
Our starting point is the single-band Hubbard Hamilto¬ 
nian: 

H =H kin + H v + H C p 
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where cj 7(T (ci ia ) is a fermionic creation (annihilation) 
operator for a particle of spin a (er =f, 4 -) at site 7 in plane 
i, rii ia = c| 70 .Ci 7cr , fV is a hopping matrix, U is the on¬ 
site Coulomb repulsion, y, is the chemical potential. We 
work with a simple cubic lattice and assume that only 
neares-neighbor hopping takes place, and that the intra- 
and inter-plane hopping parameters are equal, i.e. all 
hopping matrix elements vanish except £“ 5> = i 77 1 = 
t, where the angular brackets indicate that sites 7 and <5 
are next neighbors. We take t = 1, which sets the energy 
scale in the problem. In order to include the effect of an 
externally applied electric field, we include the following 
term: 

Hep — y ^ c7^, (2) 
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where Vi = —Vt(i/(L — 1) — 1/2) is the on-site potential 
for plane i (V is the overall potential drop across the slab 
and L is the number of planes in the slab). This form 
assumes that the charge in the slab does not screen the 
externally applied field and neglects any Coulomb inter¬ 
action between the layers themselves. We also include 
the effects of screening by including a correction to the 
bare potential of the form : 37,39 

11 Coulomb,i = kX y ' (j*;j Hbulk)\i j |> (3) 


where ribuik is the number of electrons per site in the neu¬ 
tral solid, rij = rij^ + is the average number of elec¬ 
trons at sites in plane j (site index 7 is suppressed, since 
density does not vary within a plane), and a determines 
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the interaction strength between charges, a = 2e e . In 
the last expression e is the elementary charge, e r is the 
relative dielectric constant of the solid, d is the distance 
between planes in our thin film, and A is the area per 
atom within each plane. For example, taking d = 10~ 9 m, 
A = d 2 , and e r = 10, one has a ~ 2eV. 

In the absence of screening, we use the standard DMFT 
loop: starting from a guess for the slab self-energy, 
T, a (iuj n ) = Diag(S <T i(*w n ), E cr2 (*w„),...), we obtain the 
bare Green’s function (also referred to as the Weiss 
field) of the effective action for each layer in the slab, 

Qo,itr(ikJ n )l 

^O t ia(^ n ) = \(-*Hcr(iU] r h)\ T S cr i(iu; n ), (4) 

where = (2n + 1)tt/(3 are the fermionic Matsub- 
ara frequencies at temperature T = 1//3, the index 
i = 1,2 ,... L, and Gu„ is a diagonal element of 
the fc-integrated slab Green’s function corresponding to 

Yjfj(i(jj n )'. 

GcUiOn) = —— Y -^^-, 

Nk k£BZ + m) 1 - H a (k) - E a (iu>n) 


where 1 is a L x L unit matrix (L is the number of layers 
in the slab), N k = N ka . x N ky , -/V^ and N ky are number of 
fc-points in the x and y directions of the Brillouin zone, 
respectively, and H a ( k) is a tri-diagonal L x L matrix 
in a mixed (Fourier/real-space) basis obtained from the 
sum of H k in and Hep' 
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e k = —2t(cosk x + cos fey) being the in-plane disper¬ 
sion. The impurity Green’s function for each layer, 
Gimp,icr(ikJn)t is obtained by solving the respective im¬ 
purity problem determined by Goy a . From the impurity 
Green’s function a new T, ai (iu} n ) is obtained using Eq. 
(4) with G iia replaced by G imp . icr . G a (iu n ) is only recal¬ 
culated once the solver has swept through all the layers 
and all are updated. This is repeated until conver¬ 
gence. 

To include screening, we adjust the potential according 
to Eq. (3) after each solver sweep (just before the recalcu¬ 
lation of Gimp,ia(ikJn))- For the doped calculations (not 
screened) it is also necessary to adjust the chemical po¬ 
tential of the slab between solver sweeps in order to keep 
the density fixed. 

To solve the resulting DMFT equations we use the hy¬ 
bridization expansion continuous-time quantum Monte 
Carlo solver provided as a part of TRIQS (a toolbox for 
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research on interacting quantum systems) 40_4i . In the 
segment picture 42 the on-site densities are obtained di¬ 
rectly from the solver. We estimate the quasiparticle 
weight Z from the imaginary frequency self-energy : 44 

ImSOoU " 1 
w 0 ) 

To estimate A(io = 0), the density of states at the Fermi 
level, we use the relation : 44 

/3G(/3/2) -ttA(O). ( 6 ) 

Aside from the full DMFT calculation described above, 
we also calculate the expected density distribution in the 
slab assuming that the density of electrons per site in 
each layer, n*, is determined solely by the local value of 
the chemical potential and independently of the rest of 
the layers. Within this approach, we assume that 

rii = n 3D (iii), ( 7 ) 

where n 3 D(tj) is the density per site of the three- 
dimensional Hubbard model on a cubic lattice with 
nearest-neighbor hopping at chemical potential /i and the 
same values of U, t and (3 as in the slab, and /x, = /x + v t . 
ti 3 d{h) is obtained by single-site DMFT for a set of pre¬ 
determined values of /x. The density for values of /x that 
do not coincide with values in that predetermined set is 
obtained from a linear interpolation: 

n(n) = n(nieft)x + n(n r ight)( 1 - x), (8) 

where [Meft,right are the two values of the chemical po¬ 
tential in the set which are closest to /i (such that 

l^left 33 M — fright): and X — (fright M )/(fright M/e/x)- 

In a similar manner we obtain ILA estimates for the den¬ 
sity of states at the Fermi level, A(ui = 0), and the QP 
weight Z. Two stable solutions of the DMFT equations 
of the 3D Hubbard model, a metallic and an insulating 
one, exist in parts of the U, /3, /x-range which we explore. 
For the ILA we use the density that correpsonds to the 
metallic branch. 

In all our calculations we enforce a paramagnetic 
phase, i.e. we set E-j- = Ej.. 

III. RESULTS 

A. Coexistence 

The dependence of on-site density on chemical poten¬ 
tial for the Hubbard model on a three-dimensional tight- 
binding lattice is shown in Fig. 1 for U/t = 13.2. For 
all temperatures shown there is a plateau in the n(/x) 
curve. This is due to the existence of a gap in the den¬ 
sity of states (DOS) of the solid at that value of the 
Hubbard interaction U. When the value of the chemical 


potential is inside the gap, far from the upper and lower 
Hubbard bands, the charge density does not change no¬ 
ticeably when the chemical potential is varied, due to 
the vanishing DOS in the gap. As the chemical potential 
approaches either of the two Hubbard bands, a quasipar¬ 
ticle peak appears at the Fermi level. This redistribution 
of spectral weight causes a change in the particle denisty 
and determines the end of the flat region in the n(/x) 
curve. 

As is well known, there is a region in the /i — U — T 
phase diagram of the three-dimensional fermionic Hub¬ 
bard model in which both a metallic and an insulating 
phase can exist . 45,46 When the values of /x, U, and T are 
such that the system is in this coexistence region, the 
DMFT converges to either one of the two solutions, de¬ 
pending on the initial guess (seed) for the self-energy. At 
the value of the Hubbard repulsion that we are consider¬ 
ing, U/t = 13.2, and when the system is far away from 
half-filling (when /x — [//2 is large), only the metallic so¬ 
lution exists. If one gradually brings the system closer to 
half-filling (by varying /x in small increments) and uses 
the converged value of the self-energy at each point as a 
seed for the subsequent simulation, the system remains 
on the metallic branch until the metallic solution ceases 
to exist in the region very close to half-filling. Continuing 
the recursion beyond half-filling traces out the insulating 
branch. This is reflected in the hysteresis in Fig. 1. Co¬ 
existence is possible in the slab as well. 4 ' 

In Fig. 2(a) we show the variation of the charge density 
of the metallic solution across the slab calculated with 
DMFT (symbols) and in the ILA (lines) for a number of 
different values of the applied external field at /? = 30. 
The interaction strength (U/t = 12) is chosen such that 
the bulk system at half-filling (/x = U / 2 ) would be in the 
U — T coexistence region, close to its upper boundary, 
U c 2 , beyond which only the insulating phase is stable. A 
comparison of the two approaches shows that the DMFT 
and ILA results are in good agreement in the central parts 
of the slab, whereas near the surface (within approxi¬ 
mately the first four layers) the discrepancy is significant. 
For larger fields the discrepancy between the DMFT and 
ILA is significant in a narrower region close to the surface. 
The main reason for this is the existence of highly-doped 
regions (close to the surfaces) in which the correlation 
length is smaller , 22 which limits the discrepancy between 
the two approaches to only the surface layer. We expect 
that the discrepancy between the DMFT and ILA will be 
larger whenever the correlation length is larger, in par¬ 
ticular at very low temperatures and when the first order 
nature of the transition is weaker. The same trends were 
observed earlier in calculations based on the Gutzwiller 
approximation . 21 

In Fig. 2(b) we plot the double occupancy, ( d/} = 
(n^rii]/), across the slab at zero applied external field for 
both the metallic and insulating solutions. In both the 
metallic and the insulating case the double occupancy 
is supressed at the surface. The suppression is due to 
the surface reduction of the site coordination number, 
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FIG. 1. (Color online) Bulk (3D) electron density as a 
function of chemical potential for different temperatures at 
U/t = 13.2. This value of the interaction is large enough for 
a gap to open in the spectrum of the solid. A corresponding 
plateau appears in the n(/r) curve. The simulations for dif¬ 
ferent /i values are carried out recursively. Dashed (full) lines 
correspond to down- (up-) sweep. The hysteresis indicates 
the coexistence of two solutions. Offsets added for clarity. 


which leads to an enhancement of correlation effects. The 
surface effect is also visible in the density redistribution in 
the presence of an electric field (Fig. 2(a)): the maximum 
charge deviation does not occur exactly at the surface as 
is expected from ILA approach. From Fig. 2(b) it is 
clear that in the insulating state the double occupancy 
recovers its bulk value within a single layer, whereas in 
the metallic case the recovery takes place deeper in the 
slab (approximately four layers from the surface). This 
suggests that the correlation length in the metallic state 
is longer than in the insulating state. The recovery length 
for the double occupancy at V = 0 in the metallic state 
is related to the extent of the surface effect in the case 
of applied electric field, which is given by the width of 
the discrepancy between the full DMFT calculation and 
the ILA shown in Fig. 2(a). These two lengths are of 
approximately the same magnitude. 




FIG. 2. (Color online) Charge density deviation from half- 
filling across the slab for different values of the applied field. 
Symbols indicate the full DMFT calculation; lines correspond 
to the ILA (panel a). Double occupancy across the slab for 
the insulating and metallic solution at zero field (panel b). 
U/t = 12 and /3t = 30 (both panels). 
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FIG. 3. (Color online) Electron density across the slab at 
different temperatures and different values of the electric field. 
The overall density for the slab corresponds to half-filling. 
Only "metallic" seeds were used. Symbols correspond to the 
full DMFT calculation, lines correspond to the ILA. U/t = 
13.2. 
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FIG. 4. (Color online) QP weight across the slab for the 
same simulation parameters as shown in Fig. 3. Symbols 
correspond to the full DMFT calculation, lines correspond to 
the ILA. 


B. Half-filled case 


The charge redistribution caused by an electric field 
applied to a slab at half-filling is shown in Fig. 3 
{U/t = 13.2). The value of U is outside the U — T coex¬ 
istence region for bulk (at /t = U/2 ) and is such that the 
slab is insulating at zero field. As discussed in connection 
with Fig. 1, the metallic and insulating phase can coex¬ 
ist in bulk at this value of U away from half-filling. The 
two phases can coexist in the slab as well, and we find 
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that recursive simulations for different values of V lead 
to different solutions, i.e. hysteresis (not shown). When¬ 
ever two distinct solutions exist, in Fig. 3 we show only 
the metallic one. It is clear from the figure that when 
the field applied to the slab is below a certain threshold, 
there is no redistribution of charge in the slab. When the 
field increases sufficiently, the charge density close to the 
surfaces of the film starts deviating from half-filling. The 
deviation is symmetric with respect to the center of the 
slab, as can be expected from the particle-hole symmetry 
of the Hubbard model with nearest-neighbor hopping. 

The corresponding quasiparticle weights, plotted in 
Fig. 4, are also enhanced close to the surface in the 
presence of an electric field. In contrast to what is ob¬ 
served for the density, the QP weight is enhanced close 
to the surface even at the smallest non-zero value of the 
applied field which we consider here (V = 1.5). In most 
cases, the highest value of Z in the slab occurs not imme¬ 
diately at the surface but one layer into the slab. Deeper 
into the slab Z decreases from that maximum value. At 
the highest temperature in the series {fit = 10) there is 
only a narrow central region where the enhancement is 
not significant. The width of the low-QP-weight region 
increases as the temperature is lowered. 

A{u> = 0), the density of states at the Fermi level, is 
plotted in Fig. 5 for the same simulation parameters as 
in Figs. 3 and 4. Similarly to the QP weight, A(w = 0) 
is enhanced most strongly close to the surfaces. In con¬ 
trast to what is observed for the QP weight, however, 
the enhancement of A(u = 0 ) at the lowest non-zero 
bias (V = 1.5) is insignificant. Furthermore, except at 
fit = 10, at large enough bias, the value of A(ui = 0) 
remains nearly constant (close to its maximum value in 
the slab) in the layers immediately near the surface. The 
transition to the low-DOS region in the center of the slab 
is relatively abrupt, compared with what we observe for 
the QP weight. Contrary to the width of the low-QP- 
weight region, the width of the low-DOS region decreases 
with decreasing temperature. Thus, at the highest tem¬ 
perature {fit = 10) there is a region of high Z and low 
A(ui = 0 ) close to the center of the slab, whereas at 
lower temperatures there is a region of strongly enhanced 
A(u) = 0) and weakly enhanced QP weight. 

In combination with the charge redistribution in the 
slab, the field enhancement of the QP weight and den¬ 
sity of states at the Fermi level leads to the formation 
of conductive channels close to the surfaces of the slab, 
while the central portion of the slab remains close to its 
(insulating) V = 0 state. 

The fact that a certain minimal field is necessary to 
achieve a "breakdown" of the insulator (i.e. to create an 
electrically doped region close to the surface in which the 
QP weight and DOS are significant) can be understood 
more easily by considering the dependence of density on 
chemical potential in the 3-dimensional Hubbard model 
(Fig. 1). Due to the presence of a gap in the spectrum 
for large enough U, a plateau analogous to the one in 
Fig. 1 appears in the n vs. layer curves in Fig. 3, which 


corresponds to values of the local chemical potential that 
are well within the gap. When the externally applied field 
is large enough, the local value of the chemical potential 
close to the surface of the slab deviates sufficiently from 
U /2 to induce a change in the local charge density. 

Some of our simulations (not shown) also indicate that 
when fit > 30 and the applied external field is small, only 
an insulating solution exists in the slab. In other words, 
independent of the seed, the simulation converges to the 
insulating solution throughout the slab. Above a certain 
value of the field (V ~ 3.0), two solutions reappear. This 
is consistent with what we observe in bulk (Fig. 1). 

The temperature dependence of the bulk n(ji) curves 
also elucidates what is observed in the slab. At U/t = 
13.2 the range of /Lt-values for which both a metallic and 
an insulating solution exist is clearly larger at lower tem¬ 
perature, extending to values of fi closer to the middle 
of the gap. On the metallic branch (Fig. 1), the plateau 
in nffj,) shrinks as temperature drops and almost com¬ 
pletely disappears at fi = 90. The insulating solution, 
in contrast, is less temperature sensitive. This tempera¬ 
ture dependence of the bulk metallic solution is consis¬ 
tent with what we observe in the slab (see Fig. 3): the 
minimum field strength required to "break" the insula¬ 
tor decreases as the temperature drops. The width of 
the insulating region that remains at half-filling in the 
center of the slab for a given field also decreases as the 
temperature is lowered. 

Similarly to what is observed in bulk, changes in the 
local value of the chemical potential in individual layers 
in the slab do not lead to a significant change in density, 
unless the change is large enough. 

Maximum-entropy method (MEM) 

reconstructions 48,49 of the real-frequency spectral 
function A(w) provide additional insight into the 
changes that occur across the slab as bias is increased 
and parts of the slab become metallic (see the left 
panel in Fig. 6 ). We validated our MEM results in two 
different ways: comparison with the results of Eq. ( 6 ) 
reveals very good agreement. Furthermore, calculating 
the average density in each layer from the real frequency 
spectral function, according to: 

/ +oo 

dkjAi(u))f(u - Hi), (9) 

-OO 

where A,(cc) is the MEM spectral function of layer i and 
f(uj) is the Fermi-Dirac distribution function, we find 
excellent agreement with the Monte-Carlo charge density 
data (Fig. 1). The appearance and development of the 
quasiparticle peaks due to the spectral transfer from the 
nearby Hubbard bands (the frequency-integrated DOS 
remains constant) is clearly visible in the figure. The 
width of the peaks is proportional to the QP weight Z 
and reaches a maximum approximately one layer into 
the slab from each of the surfaces, which is consistent 
with Fig. 4. Similarly, the abrupt rise of A(ui = 0) is 
reproduced (in the MEM plot A{u> = 0) corresponds to 
the height of the QP peak). The MEM results for the 
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breakdown bias. The differences between the ILA esti¬ 
mates for Z and A(0) and the full calculation are more 
pronounced. The largest discrepancy in Z is observed at 

V = 3, the bias value at which the system is closest to 
the breakdown point. At V = 3 and fit > 30 the ILA 
and full DMFT results for Z differ not only near the sur¬ 
face, but through a significant part of the slab. A similar 
observation is valid for A(uj = 0), but in addition to the 
significant disagreement between the ILA and DMFT es¬ 
timates of A(u = 0) at V = 3, the ILA also significantly 
overestimates the value of d(w = 0) at V = 1.5 and 
fit = 60, 90. The agreement is better at larger values of 

V and at higher temperature (lower /3). 


C. Screening 


FIG. 5. (Color online) The variation of the DOS at the Fermi 
level across the slab corresponding to the density and QP 
weight data shown in Figs. 3 and 4. Symbols correspond to 
the full DMFT calculation, lines corrspond to the ILA. 


doped case (Fig. 6, right) are discussed in more detail in 
the following section. 
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FIG. 6. (Color online) Maximum-entropy method reconstruc¬ 
tion of the spectral functions Ai(u}) throughout the slab. 
Data is interpolated between layers for clarity. Left panel: 
U/t = 13.2,/x = U/2, fit = 30, V = 4.5 and average density 
corresponding to half filling. Right panel: U/t = 16.2, fit = 
10, V = 6 and average density per site in the slab n = 1.04. 
The MEM results are in very good agreement with those ob¬ 
tained by Eq. (6). 


We also investigate the effect of screening on the half- 
filled slab (see Fig. 7). In Figs. 7(a,c) we show the 
density and potential across a slab in the metallic phase 
(U/t = 10 and fit = 10) in the presence of screening 
and in the unscreened case. For the screened case we 
use a = 1. The total potential drop across the slab in 
both cases is the same (this means that in the screened 
case the applied external field V is greater). In prac¬ 
tice, we first fix the externally applied field V and the 
screening strength a and allow the simulation to con¬ 
verge to a self-consistent solution that satisfies Eq. (3). 
This leads to a reduced potential drop across the slab. 
Then a non-screened simulation with the same potential 
drop is performed for comparison. The effect of screening 
is to reduce the penetration of the field in the slab. As 
a consequence, the density deviation (from half-filling) 
is reduced, especially in the central region of the slab. 
Qualitatively the results remain very similar: the devi¬ 
ation from half-filling is symmetric with respect to the 
central region and largest close to the two surfaces in the 
slab. Screening masks the influence of the surface on the 
density: in the absence of screening the maximum of the 
denisty does not occur at the surface but one layer into 
the slab. In contrast, in the screened case, the maximum 
is at the surface. The surface reduction effect is still visi¬ 
ble, but softened, due to the decreased penetration of the 
field in the slab. The effect of screening on an insulating 
slab ( U/t = 13.2, fit = 30) is analogous, Figs. 7(b,d). In 
contrast to the metallic case, however, the effect is much 
less pronounced. 


It is clear (Fig. 3) that the ILA estimate of the density 
agrees fairly well with the full calculation for almost all 
the temperatures and field strengths we consider. The 
difference between the two approaches is greatest close 
to the surfaces of the slab, where the independent plane 
approximation tends to overestimate the density and fails 
to capture the dip evident in the full calculation, and for 
values of V which are close to the "breakdown" bias of 
the insulator. The ILA underestimates the value of the 


D. Doped case 

The effect of the field on a doped slab is quite different. 
In Figs. 8 - 10 we show the charge density redistribution 
(DMFT and ILA results), A(ui = 0), and quasiparticle 
weight across the slab for four different values of U and 
f) ( U/t = 13.2,16.2; fit = 10,30). At our chosen doping 
(n = 1.04) and in the absence of a field (V = 0), both the 
density of states at the Fermi level (Fig. 9) and the QP 
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FIG. 7. (Color online) The effect of screening on a metallic 
slab, panels a, c, and an insulating slab, panels b, d. The 
effect is much more pronounced for the metallic case, a is the 
screening parameter described in Sec. II. 


weight (Fig. 10) are significant throughout the slab in all 
four cases, indicating a metallic phase. In this case even a 
small field causes charge redistribution. Since the average 
density in the slab is kept fixed, the pile-up of charge on 
one side is accompanied by a decrease in the charge den¬ 
sity on the other side. However, if the value of U is large 
enough (i.e. sufficient for a gap to appear in the DOS 
of the solid at half-filling), this decrease in the charge 
density on the depleted side pauses as soon as half-filling 
is reached. From that point on, as the magnitude of V 
increases the charge build-up on the other side is com¬ 
pensated by an increase in the width of the half-filling 
region. This is what we observe here for U/t = 13.2 and 
U/t = 16.2. Thus, at large U, as the field increases one 
of the sides of the slab becomes more conducting while 
an insulating layer of increasing thickness develops at the 
opposite end of the slab. Note that the density redistri¬ 
bution is accompanied by a change in the QP weight: the 
QP weight is enhanced on the side with excess charge and 
suppressed on the "depleted" side. The density of states 
at the Fermi level follows a similar trend. As the temper¬ 
ature is lowered from fit = 10 to fit = 30 the transition 
between the low- and high-conductivity regions becomes 
more abrupt (Fig. 9). Variation of U/t in the range con¬ 
sidered (13.2 — 16.2) does not significantly influence the 
density results. The effect of increasing U on A(u> = 0) 
is most pronounced at V = 0, in which case it leads to 
a decrease of the DOS at the Fermi level throughout the 
slab. 

The insets in Fig. 8 show how the chemical potential 
should be adjusted in order to keep the average amount 
of charge per layer in the slab constant (at n = 1.04) at 
each bias. The adjustment is necessary due to the non¬ 
linear dependence of charge density on /i (see Fig. 1). 
Experimentally, this situation corresponds to a slab that 
is electrically insulated from its environment (the overall 
amount of charge in the slab remains constant). The 
observed reduction of chemical potential in the presence 
of external electric field indicates that if connected to 
a charge reservoir (e.g. via metallic leads) an overdoped 


system would tend to absorb charge when an electric field 
is applied. In simulations with doped slabs at constant 
lx we do observe an increase in the amount of charge in 
the slab as the field increases (not shown). This charge 
absorption is cuased by the nonlinearity of the n(/r) curve 
and not by a potential difference between the system and 
its environment. 
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FIG. 8. (Color online) Variation of the electron density 
throughout a 24 layer slab for selected values of the electric 
field, Hubbard U and temperature. The average density for 
the whole slab is fixed at 1.04 electrons per site. The insets 
show how the chemical potential of the slab changes as the 
field is increased. This adjustment is necessary to keep the 
average electron density in the slab fixed. 
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FIG. 9. (Color online) Variation of the density of states at the 
Fermi level across a 24 layer slab for the same values of the 
electric field, Hubbard U and temperature as in Fig. 8. The 
average density for the whole slab is fixed at 1.04 electrons 
per site. Symbols correspond to the full DMFT calculation, 
lines corrspond to the ILA. 




























































At higher values of the interaction {U/t = 16.2) oscil¬ 
lations appear in the QP weight vs. layer number. The 
effect of lowering temperature is to make the oscillations 
more pronounced (cf. Fig. 10, lower two panels). In the 
presence of a significant DOS, these oscillations would 
correspond to alternating regions of low and high mo¬ 
bility of the quasiparticles in the film. However, since 
A(u) = 0) is rather low in the part of the film where these 
oscillations occur (see Fig. 9), in practice the variation 
in conductivity would be small. The fact that these os¬ 
cillations in Z are also reproduced by the ILA indicates 
that they are not related to the slab geometry. 

It is worth noting that different parts of the slab can 
coexist in different states, i.e. that the metallicity of 
one part of the slab does not penetrate throughout the 
slab to destroy the low-QP-weight and low-DOS region 
in other parts of the slab. The MEM reconstruction of 
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FIG. 10. (Color online) Variation of the QP weight across a 
24 layer slab for the same values of the electric field, Hubbard 
U and temperature as in Figs. 8 and 9. The average density 
for the whole slab is fixed at 1.04 electrons per site. Symbols 
indicate the full DMFT calculation; lines correspond to the 
ILA. 

the slab spectral functions is shown in the right panel 
of Fig. 6 for U/t = 16.2, fit = 10, and V = 6. In this 
case the agreement with Eq. (6) and DMFT density data 
is again very good. Due to the doping, the position of 
the lower and upper Hubbard bands with respect to the 
Fermi level is not symmetric, in contrast to the half-filled 
case. This asymmetry is also reflected in Fig. 8 to Fig. 
10. From the MEM reconstruction it is clear that once 
again the QP and DOS enhancement on the right side of 
Figs. 10 and 9 are associated with the appearance of a 
QP peak. On the basis of the MEM reconstruction, it is 
to be expected that at sufficiently strong fields the lower 
Hubbard band will also approach the Fermi level and a 
quasiparticle peak will appear on the low-doped side as 
well, which will cause a "secondary breakdown" to occur 
at the insulating side of the film. 


The doping and the resulting asymmetry lead to a very 
different response compared to the half-filled case. In the 
latter case, the application of the field to the insulating 
slab leads to the formation of symmetric zones of rela¬ 
tively high conductivity close to the surfaces, whereas in 
the doped case, which is metallic in the absence of a field, 
application of a field causes one side to become more con¬ 
ductive, whereas the other side becomes insulating, until 
the field is large enough for the "secondary" breakdown 
to occur. 

The ILA results for the charge density match the full 
calculation very nearly in almost all cases, except in the 
surface layers, and - when the bias is in the vicinity of 
the "secondary breakdown" - in a narrow region close 
to the surface on the insulating side. The agreement 
on the over doped side improves as the bias increases. 
We attribute this to the short correlation length in the 
highly-doped (metallic) regime. On the underdoped side 
the discrepancy is most significant at low temperature 
{fit > 30), close to the "secondary breakdown". As in 
the half-filled case, the ILA underestimates the break¬ 
down voltage. This tendency to overestimate the extent 
of the metallic phase is reflected also in the ILA results 
for the density of states, Fig. 9, and QP weight, Fig. 10. 
Overall, the agreement between the DMFT and ILA is 
much better for the doped case than for the half-filled 
slab, which can be partially attributed to the fact that 
at most biases and temperatures a large part of the slab 
is metallic. 

When the Hubbard interaction is lowered to U/t = 
10.2 (Fig. 11), no half-filling region of significant width 
develops, as there is no plateau in the n(/i) curve at this 
interaction strength. At this value of U the QP weight 
exhibits a minimum as a function of layer number, and 
is enhanced close to the slab surfaces. Expectedly, as in 
the large U case, even a weak field is sufficient to cause 
a redistribution of charge. The agreement between the 
DMFT and the ILA results for density, A(ui = 0), and 
QP weight is excellent at all bias values we consider. The 
largest discrepancy between the two approaches occurs 
in the two layers immediately at the surface and is most 
significant in the value of Z (panel c). 


IV. CONCLUSION 

In summary, we investigate the properties of strongly 
correlated thin films under bias using IDMFT and ex¬ 
amine the validity of a computationally cheaper approx¬ 
imation. We observe switching behavior in both half- 
filled and doped films. In the half-filled slab a sufficiently 
strong field (larger than some threshold) is necessary to 
produce conducting regions near the surfaces of the film. 
For doped films, there is no threshold field and the ap¬ 
plication of a field initially causes one side of the film to 
become more insulating, before a secondary breakdown 
occurs. Taking screening into account does not lead to 
qualitative changes of the results. We therefore expect 
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FIG. 11. (Color online) Electron density, DOS at the Fermi 
level, and QP weight across a doped 24 layer slab for the same 
values of the electric field as in Figs. 8 - 10, but lower U. The 
average density in the slab is fixed at 1.04 electrons per site. 
The inset shows the change of fi necessary to keep the average 
density constant for different values of V. Symbols correspond 
to the full DMFT calculation, lines corrspond to the ILA. 


that the conclusions for the non-screened case will re¬ 
main valid in the presence of screening. The hysteretic 
behavior associated with the first order Mott transition 


observed in bulk persists in the slab, and should lead to 
memory effects in devices. 

In spite of the breakdown of useful concepts from band 
theory, such as band bending, the local independent layer 
approximation accurately reproduces the full IDMFT 
calculations in both the half-filled and doped cases, ex¬ 
cept in the layers immediately at the surface of the slab 
and close to transition points. Our calculations confirm 
and extend earlier findings in theoretical studies of colos¬ 
sal magnetoresistance, where DMRG calculations in ID 
systems indicated the existence of a universal density- 
potential relation of the interface Mott transition. 50 Thus 
the ILA may be used as a fairly reliable and quick first 
estimate for the charge density (and to a lesser extent 
for the quasiparticle weight and DOS at the Fermi level). 
This may allow for calculations of device properties such 
as e.g. electrostatic charge distribution, switching be¬ 
havior, differential capacitance, which are not easily ac¬ 
cessible with full DMFT calculations due to the larger 
computational cost. 
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